function [x k] = gauss_seidel(a,b, err_max, iter_max)
% Gauss-Seidel method for solving for x when b=ax
% This method will only converge if the matrix a is either:
% 1. diagonally dominant (SDD) or
% 2. symmetric and positive definite
% Argurments
% a: square MxM matrix
% x,b: Mx1 column vectors
% err_max: maximum error tolerance
% iter_max: maximum number of iterations until stopping

if nargin < 3
  err_max = 10^-7;
  iter_max = 15;
elseif nargin < 4
  iter_max = 15;
end

[m n] = size(a);
[br bc] = size(b);

if m ~= n
  msg = sprintf('ERROR:\tmatrix a must be square %d x %d', m,m);
  disp(msg);
  x = NaN;
  return
elseif br ~= m
  if bc == m & br == 1
    msg = sprintf('WARNING:taking transpose of b to transform it to %d x 1', m);
    disp(msg);    
    b = b' 
  else
    msg = sprintf('ERROR:\tb is wrong size (%d x %d) should be %d x %d', br, bc, m, 1);
    disp(msg);        
    x = NaN;
    return;
  end
end
    

l = tril(a,-1);    % lower triangular matrix 
u = triu(a,1);     % upper triangular matrix
d = diag(diag(a)); % diagonal matrix

% inverses of triangular matrices are trivial
ld_inv = inv(l+d);
ld_inv_b = ld_inv * b;
ld_inv_u = ld_inv * u;

x= zeros(m,1);     % initialize iterative solution
  
for k=1:iter_max
  x_next = ld_inv_b - ld_inv_u*x;
%  for i=1:m
%    x(i) =  (1/a(i,i))*(b(i) - sum(a(i,i+1:m) .* x_old(i+1:m)) - ...
%                        sum(a(i,1:i-1) .* x_old(
%  end
  if abs(norm(x-x_next)) < err_max
    break; 
  end
  x = x_next;  
end
